;;----------------------------------------------------------------------------- 
;; ncl script for post processing 
;; written by Kai Zhang, MPI-M, 2009 
;; name: mkdata_mm.ncl 
;;----------------------------------------------------------------------------- 

;;-----------------------------------------------
;; load libraries  
;;-----------------------------------------------

  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
 
;;-----------------------------------------------
;; set model name  
;;-----------------------------------------------

  load "tmp_env.ncl" 
  load "tmp_rcase.ncl" 

  L_CAM = 1
  L_ECHAM = 1
  L_GAMIL = 1

  if (vmod.eq."cam")
     L_CAM = 0
  end if

  if (vmod.eq."echam")
     L_ECHAM = 0
  end if

  if (vmod.eq."gamil")
     L_GAMIL = 0
  end if

  L_debug = 1


;;-----------------------------------------------
;; file name    
;;-----------------------------------------------
  
   inname = "out_"+vtype+"_"+vmod+".nc"
   innama = "out_PS_"       +vmod+".nc" 
   innama = "out_met2d_"    +vmod+".nc" 
   innamb = "out_T_"        +vmod+".nc" 
   innamb = "out_met3d_"    +vmod+".nc" 
   innamc = "out_postmm_"   +vmod+".nc"

   outnama = rcase + "." + varname + "Z." + year + "." + vmod + ".nc"
   system("rm " + outnama )
   fileoua = addfile(outnama,"c")
 
   outnamb = rcase + "." + varname + "P." + year + "." + vmod + ".nc"
   system("rm " + outnamb )
   fileoub = addfile(outnamb,"c")
 
   outnamc = rcase + "." + varname + "S." + year + "." + vmod + ".nc"
   system("rm " + outnamc )
   fileouc = addfile(outnamc,"c")

   filein = addfile(inname,"r")
   fileia = addfile(innama,"r")
   fileib = addfile(innamb,"r")
   fileic = addfile(innamc,"r")

   print(""+inname)
   print(""+innama)
   print(""+innamb)
   print(""+innamc)

;;-----------------------------------------------
;; file dims 
;;-----------------------------------------------

   if (L_ECHAM.eq.0) then
      T = fileic->st
   end if   

  dims = dimsizes(T) 

  ntim = dims(0)
  nlev = dims(1)
  nlat = dims(2)
  nlon = dims(3)   

;;-----------------------------------------------
;; file attribute  
;;-----------------------------------------------

   filedimdef(fileoua,"time",-1,True)
   filedimdef(fileoub,"time",-1,True)
   filedimdef(fileouc,"time",-1,True)

;;-----------------------------------------------
;; cycle of vars  
;;-----------------------------------------------

   do ia = 0, dimsizes(varname) - 1
 
   
   if (L_ECHAM.eq.0) then
      vart = filein->$vecname(ia)$
      surp = fileia->aps
      T = fileic->st
   else
      vart = filein->$varname(ia)$
      surp = fileia->PS
      T = fileib->T
   end if
 


 
   L_var_interface = 1
 
   if (varname.eq."CMFMC") then
      L_var_interface = 0
   end if 
 
   if (L_var_interface.eq.0) then
      levp = dimsizes(vart(0,:,0,0)) 
      varx = vart(:,0:levp-2,:,:) 
   else
      varx = vart
   end if 


   if (varname.eq."CLOUD") then
      varx@units = "percent" 
   end if


 
  I_lev_14 = 0 
  I_lev_10 = 1 

  if(I_lev_14.eq.0) then 
     varp = new( (/ ntim, 14, nlat, nlon /), float )
  end if 

  if(I_lev_10.eq.0) then 
     varp = new( (/ ntim, 10, nlat, nlon /), float )
  end if



  
  pmtop = 2.194 
  
  pold  = new( (/nlev,nlat,nlon/), float ) ;hpa


;;-----------------------------------------------
;; load vertical coordinate 
;;-----------------------------------------------

  load "mkdata_mm_vpara.ncl"


  do itim = 0, ntim - 1  


  ;; gamil model 


  if (L_GAMIL.eq.0) then

  if (L_var_interface.eq.0) then

    do k = 0, nlev - 1 
      pold(k,:,:) = ( surp(itim,:,:)/100. - pmtop ) * sigmi(k)/1000. + pmtop    ; hpa
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do

  else

    do k = 0, nlev - 1 
      pold(k,:,:) = ( surp(itim,:,:)/100. - pmtop ) * sigma(k)/1000. + pmtop    ; hpa
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do

  end if 

  end if


  ;; cam model 

  if (L_CAM.eq.0) then

  if (L_var_interface.eq.0) then

    do k = 0, nlev - 1 
      pold(k,:,:) = ( surp(itim,:,:)/100.) * (hybi_cam(k)) + 1013.25 * (hyai_cam(k))
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do

  else

    do k = 0, nlev - 1 
      pold(k,:,:) = ( surp(itim,:,:)/100.) * (hybm_cam(k)) + 1013.25 * (hyam_cam(k))
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do

  end if

  end if

  ;; echam model 

  if (L_ECHAM.eq.0) then

  if (L_var_interface.eq.0) then

    do k = 0, nlev - 1 
      pold(k,:,:) = ( surp(itim,:,:)/100.) * (hybi_echam(k)) + (hyai_echam(k)) / 100.
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do

  else

    do k = 0, nlev - 1 
      pold(k,:,:) = ( surp(itim,:,:)/100.) * (hybm_echam(k) ) + ( hyam_echam(k) ) / 100.
      if (L_debug.eq."0")
          print("MAX,MIN,AVG of Pressure = "+max(pold(k,:,:))+"   "+min(pold(k,:,:))+"   "+avg(pold(k,:,:)))
      end if
    end do

  end if

  end if
 


;...................................................................................
; unit
;...................................................................................

  R = 8.3143

  if (vtypa.eq."num") then
     varx(itim,:,:,:) = varx(itim,:,:,:) * 1013.25 * 100. * 29. / ( R * 273.15 * 1000. )    ; kg / m3
     varx(itim,:,:,:) = varx(itim,:,:,:) * 1.E-06                   ; number / cm3
     varx@units    = "N cm~S~-3~E~"
  end if

  if (vtypa.eq."rd".or.vtypa.eq."rw") then
     varx(itim,:,:,:) = varx(itim,:,:,:) * 1.e+06
     varx@units = "um" 
  end if

  varx@long_name = varname(ia)


 ; do k=0,nlev-1 
 ;   pres(k,:,:) = surp*sigma(k)/1000./100.
 ; end do

 ;print(pres)



  I_SO4_molecules_cm3 = 0 

  ltmp = varname.eq."SO4"      .or. \
         varname.eq."SO4_gas" 

  lyes = ltmp .and. I_SO4_molecules_cm3.eq.0  

  if (lyes) then 

     if (L_ECHAM.eq.0) then
        varx(itim,:,:,:) = varx(itim,:,:,:) * 3. ;kg(S)->kg 
     end if 

     ;;do ik = 0, dimsizes(varx(0,:,0,0)) - 1
     ;;   varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 1e12 * 29. / 96.  ;kg(SO4)/kg -> ppt
     ;;end do

     ;;varx@units = "pptv"

     do ik = 0, dimsizes(varx(0,:,0,0)) - 1  
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 1000. / 96.                                                       ;kg/kg -> mol(SO4)/kg(air)
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * pold(ik,:,:) * 100. * 29. / ( 8.314409 * T(itim,ik,:,:) * 1000.)     ;mol(SO4)/m3(air) 
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 1.E-06 * 6.02217E+23                                              ;molecules(SO4)/cm3(air)
     end do

     varx@units = "molecules cm-3"

    ;comments
    ;in echam5-ham, the unit of reactive gases is volumn mixing ratio ( 29 * kg(SO4) / molarweight(so4) ), or ppt*1E+12 

     print("max,min of SO4 " + max(varx(itim,:,:,:)) + " " + min(varx(itim,:,:,:)) + " molecules cm-3 " ) 
     
  end if


    ;do not calculate ppt from density
    ;density of dms     0.840 g/cm3 (not vaport phase!) 
    ;                   2.100 * 1.E-03 g/cm3
    ;density of so2     2.551 * 1.E-03 g/cm3 (vapor) 


  I_DMS_pptv = 0

  ltmp = varname.eq."DMS" 

  lyes = ltmp .and. I_DMS_pptv.eq.0 

  if (lyes) then

     if (L_ECHAM.eq.0) then
        varx(itim,:,:,:) = varx(itim,:,:,:) * 62./32. ;kg(S)->kg 
     end if

     do ik = 0, dimsizes(varx(0,:,0,0)) - 1
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 1e12 * 29. / 62.  ;kg(DMS)/kg -> ppt
     end do

     varx@units = "pptv"
     
     print("max,min of DMS " + max(varx(itim,:,:,:)) + " " + min(varx(itim,:,:,:)) + " ppt " ) 

  end if



  doit = varname.eq."MSO" .or. varname.eq."MSU".or. varname.eq."MBC".or.varname.eq."MPM".or.varname.eq."MSS" .or. \
         varname.eq."MDU".or.varname.eq."MSOANT".or.varname.eq."MSOBIO" .or. \
         varname.eq."OC_KS" .or. varname.eq."OC_AS" .or. varname.eq."OC_CS" .or. varname.eq."OC_KI" 

  if ( doit ) then

     if (L_ECHAM.ne.0.and.varname.eq."MSU") then
        varx(itim,:,:,:) = varx(itim,:,:,:) / 3. ;kg(S)->kg 
        print("varx(itim,:,:,:) = varx(itim,:,:,:) / 3. for MSU") 
     end if

     do ik = 0, dimsizes(varx(0,:,0,0)) - 1

 ;unit conversion based on ideal air state equation
 ;http://www.grc.nasa.gov/WWW/K-12/airplane/eqstat.html
 ;
 ;P (Pa) V (m3) T (K) M (kg) 
 ;P * V / T = M * 1000 / 29 * 8.31 
 ;rou = M / V = P * 29 / ( T * 1000. * 8.31 ) 
;;printVarSummary(varx) 
;;printVarSummary(pold) 
;;printVarSummary(T) 
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * pold(ik,:,:) * 100. * 29 \
                             / ( 8.314409 * T(itim,ik,:,:) * 1000. )      ;kg(S)/kg(air) -> kg(S)/m3(air) 
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 1.E+09                   ;kg(S)/m3(air) -> ug(S)/m3(air)

     end do

     if (varname.eq."MSU") then
     varx@units = "ug(S)/m3"
     else
     varx@units = "ug/m3"
     end if

  end if








  I_SO2_pptv = 0

  ltmp = varname.eq."SO2" 

  lyes = ltmp .and. I_SO2_pptv.eq.0 

  if (lyes) then

  if (L_ECHAM.eq.0) then
     varx(itim,:,:,:) = varx(itim,:,:,:) * 2. ;kg(S)->kg(SO2)
  end if

     do ik = 0, dimsizes(varx(0,:,0,0)) - 1  
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 1e12 * 29. / 64  ;kg(SO2)/kg -> ppt
     end do

     varx@units = "pptv"

     print("max,min of SO2 " + max(varx(itim,:,:,:)) + " " + min(varx(itim,:,:,:)) + " ppt " ) 
     
  end if


  I_RN222_mBqm3 = 0

  ltmp = varname.eq."RN222"

  lyes = ltmp .and. I_RN222_mBqm3.eq.0

  if (lyes) then

     do ik = 0, dimsizes(varx(0,:,0,0)) - 1
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 29. / 222. * 1.e+21  
        varx@units = "mol/mol"
        varx(itim,ik,:,:)  = varx(itim,ik,:,:) * 5.6 * 1e-02 * 1.e+03  
        varx@units = "mBq/m3"
     end do


     print("max,min of RN222 " + max(varx(itim,:,:,:)) + " " + min(varx(itim,:,:,:)) + " mBq/m3 " )

  end if




  if(I_lev_14.eq.0) then 
     lev = new(14,"float")
     lev = (/10,20,50,70,100,150,200,300,400,500,600,700,850,1000/)
  end if 

  if(I_lev_10.eq.0) then 
     lev = new(10,"float")
     lev = (/100,200,300,400,500,600,700,800,900,1000/)
  end if     



  L_vinth2p = True
  

 ;print(surp)

  if (L_vinth2p) then


     if (L_ECHAM.eq.0) then
        ;;fake = hyam_echam 
        ;;tmpp = vinth2p(varx(itim,:,:,:),fake,hybm_echam,lev*100.,surp(itim,:,:)*100.,1,1.,1,False)
        fake = hyam_echam / 101325.   
        tmpp = vinth2p(varx(itim,:,:,:),fake,hybm_echam,lev,surp(itim,:,:),1,1013.25,1,False)
     end if
     
     if (L_CAM.eq.0) then
        fake = hyam_cam
        tmpp = vinth2p(varx(itim,:,:,:),fake,hybm_cam,lev,surp(itim,:,:),1,1000.,1,False)
     end if
     
     if (L_GAMIL.eq.0) then
        fake = sigma
        fake = 0.
        tmpp = vinth2p(varx(itim,:,:,:),fake,sigma/1000.,lev,surp(itim,:,:),1,1000.,1,False)
     end if
     
     if ( (L_GAMIL+L_CAM+L_ECHAM) .ne. 2 ) then
        print(" CHECK L_GAMIL, L_CAM, and L_ECHAM ! ") 
     end if
     
     print("max,min varx: "+max(varx(itim,:,:,:)) +" "+min(varx(itim,:,:,:)))    
     print("max,min tmpp: "+max(tmpp) +" "+min(tmpp))    
     
     
     varp(itim,:,:,:) = (/ tmpp(:,:,:) /)

    ;printVarSummary(varp)
    ;printVarSummary(tmpp)
  
  else
 
    do   j = 0, nlat - 1
      do i = 0, nlon - 1
        varp(itim,:,j,i) = linint1( pold(:,j,i), varx(itim,:,j,i), False, lev, 0 )
      end do 
    end do 
  	
  end if
  	 

 ;printVarSummary(varx)
 ;printVarSummary(varp)
 
 
 ;print(""+varp(itim,9,20,20))

;...................................................................................
;  
;...................................................................................
 
 do i=0,5
    varp(itim,i,:,:) = mask(varp(itim,i,:,:),(lev(i).gt.surp(itim,:,:)/100.),False)
 end do


;...................................................................................
;  
;...................................................................................




 ;delete(pres)
 ;delete(lev)


  end do 

  delete(pold)




  varp!0              = varx!0
  varp&time           = varx&time
  varp&time@units     = varx&time@units
  varp&time@long_name = "time"

  varp!1              = "lev"
  varp&lev            = lev
  varp&lev@units      = "hPa"
  varp&lev@long_name  = "Pressure (hPa)"

  varp!2              = varx!2
  varp&lat            = varx&lat
  varp&lat@units      = varx&lat@units
  varp&lat@long_name  = varx&lat@long_name

  varp!3              = varx!3
  varp&lon            = varx&lon
  varp&lon@units      = varx&lon@units
  varp&lon@long_name  = varx&lon@long_name

  varp@units          = varx@units
  varp@long_name      = varx@long_name
  varp@model          = vmod
  varp@season         = vseason
  varp@varname        = varname
  varp@type           = vtype
  varp@plot           = "pressure"
  varp@owner          = "Kai Zhang"




   shape = dimsizes(varp)
   nlev  = shape(1)
   nlat  = shape(2)

   varz =  new( (/ntim,nlev,nlat,1/), float )




; varz = dim_avg(varp)

; varx&time = 0.

  varz!0              = "time"
  varz&time           = varx&time
; varz&time@units     = varx&time@units
; varz&time@long_name = varx&time@long_name

  varz!1              = "lev"
  varz&lev            = lev
  varz&lev@units      = "hPa"
  varz&lev@long_name  = "Pressure (hPa)"

  varz!2              = "lat"
  varz&lat            = varx&lat
; varz&lat@units      = varx&lat@units
; varz&lat@long_name  = varx&lat@long_name

  varz!3 = "lon"
  varz&lon = varx&lon(0)

  varz@long_name      = "zonal mean "+varp@long_name
  varz@model          = vmod
  varz@season         = vseason
  varz@varname        = varname
  varz@type           = vtype
  varz@plot           = "zonal"
  varz@owner          = "Kai Zhang"

  varz(:,:,:,0) = dim_avg(varp(:,:,:,:))




  vars =  new( (/ntim,nlat,nlon/), float )


  vars!0              = varx!0
  vars&time           = varx&time
  vars&time@units     = varx&time@units
  vars&time@long_name = "time"


  vars!1              = varx!2
  vars&lat            = varx&lat
  vars&lat@units      = varx&lat@units
  vars&lat@long_name  = varx&lat@long_name

  vars!2              = varx!3
  vars&lon            = varx&lon
  vars&lon@units      = varx&lon@units
  vars&lon@long_name  = varx&lon@long_name


  vars@units          = varx@units
  vars@long_name      = varx@long_name
  vars@model          = vmod
  vars@season         = vseason
  vars@varname        = varname
  vars@type           = vtype
  vars@plot           = "surface"
  vars@owner          = "Kai Zhang"

  vars(:,:,:) = varx(:,dimsizes(varx(0,:,0,0))-1,:,:)


  fileoua->var=varz
  fileoub->var=varp
  fileouc->var=vars


  end do

